function [rej_rate,p_adj] = rej_rate_IM_mean(SigmaHat,membership,alpha_sig,alt)

rng(7)
n = size(SigmaHat,1);
clusters = unique(membership);
G = numel(clusters);

nSim = 10000;  % number of simulation replications

% only generates group mean
L = repmat(membership',G,1)==repmat(clusters,1,n);
N_g_vec = sum(L,2);
L_weighted = L./repmat(N_g_vec,1,n);
b_mat = alt+(randn(nSim,G)*chol(L_weighted*SigmaHat*L_weighted'))';
t_vec = nanmean(b_mat)./(nanstd(b_mat)/sqrt(G-1));
test_vec = abs(t_vec)>tinv(1-alpha_sig/2,G-1);
rej_rate = mean(test_vec);

% adjusted critical values and corresponding p-value threshold such that 
% rejection rate is alpha_sig
cv_adj = quantile(abs(t_vec),1-alpha_sig); 
p_adj = 2*(1-tcdf(cv_adj,G-1));


